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Abstract 

Common efficient schemes for the incompressible Navier-Stokes equations, 
such as projection or fractional step methods, have limited temporal accu- 
racy as a result of matrix splitting errors, or introduce errors near the domain 
boundaries (which destroy uniform convergence to the solution). In this pa- 
per we recast the incompressible (constant density) Navier-Stokes equations 
(with the velocity prescribed at the boundary) as an equivalent system, for 
the primary variables velocity and pressure. We do this in the usual way 
away from the boundaries, by replacing the incompressibility condition on 
the velocity by a Poisson equation for the pressure. The key difference from 
the usual approaches occurs at the boundaries, where we use boundary con- 
ditions that unequivocally allow the pressure to be recovered from knowledge 
of the velocity at any fixed time. This avoids the common difficulty of an, 
apparently, over-determined Poisson problem. Since in this alternative for- 
mulation the pressure can be accurately and efficiently recovered from the 
velocity, the recast equations are ideal for numerical marching methods. The 
new system can be discretized using a variety of methods, in principle to any 
desired order of accuracy. In this work we illustrate the approach with a 
2-D second order finite difference scheme on a Cartesian grid, and devise an 
algorithm to solve the equations on domains with curved (non-conforming) 
boundaries, including with a non-trivial topology (a circular obstruc- 

tion inside the domain). This algorithm achieves second order accuracy in 
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the L°° norm, for both the velocity and the pressure. The scheme has a 
natural extension to 3-D. 
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1. Introduction. 

A critical issue in the numerical solution of the incompressible Navier- 
Stokes equations is the question of how to implement the incompressibility 
constraint. Equivalently, how to recover the pressure from the flow velocity, 
given the fact that the equations do not provide any boundary condition 
for the pressure. This has been an area of intense research, ever since the 



pioneering MAC scheme [Ij] of Harlow and Welch in 1965. Of course, one 
can avoid the problem by simultaneously discretizing the momentum and the 
divergence free equations, as in the difference scheme proposed by Krzywicki 



and Ladyzhenskaya [22[, which can be shown to converge — while avoiding 
the need for any pressure boundary conditions. Approaches such as these, 
however, do not lead to efficient schemes. 

Generally the dilemma has been that of a trade-off between efficiency, 
and accuracy of the computed solution near the boundary. However, many 
applications require both efficiency, and accuracy. For example, to calculate 
fluid solid interactions, both the pressure and gradients of the velocity are 
needed at the solid walls, as they appear in the components of the stress 
tensor. Furthermore, these objectives must be achievable for "arbitrary" 
geometries, not just simple ones with symmetries that can be exploited. 
Unfortunately, these requirements are not something that current algorithms 
are generally well suited for, as the brief review below is intended to showll] 
However, we believe that algorithms based on a Presure Poisson Equation 
(PPE) reformulation of the Navier-Stokes equations — reviewed towards 
the end of this introduction — offer a path out of the dilemma. The work 
presented in this paper is, we hope, a contribution in this direction. 

Projection methods are very popular in practice because they are efficient. 
They achieve this efficiency by: (i) Interpreting the pressure as effecting a 
projection of the flow velocity evolution into the set of incompressible flelds. 



^This is not intended as a thorough review of the field, and we apologize for the many 
omissions. 
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That is, write the equations in the form = P (|U A u — (V ■ u) u + f), 
where V is the appropriate projection operator, /z is the kinematic viscosity, 
u is the flow velocity vector, and f is the vector of applied body forces, (ii) 
Directly evolving the flow velocity. The question is then how to compute V. 

In their original formulation by Chorin |6] and Temam 38|, the projec- 
tion method was formulated as a time splitting scheme in which: First an 
intermediate velocity is computed, ignoring incompressibility. Second, this 
velocity is projected onto the space of incompressible vector fields — by solv- 
ing a Poisson equation for pressure. Unfortunately this process introduces 
numerical boundary layers into the solution, which can be ameliorated (but 
not completely suppressed) for simple geometries — e.g. ones for which a 
staggered grid approach can be implemented. 

The development of second order projection methods ji], [isl, 20, 27, 39 
provided greater control over the numerical boundary layers. These are the 
most popular schemes used in practice. However, particularly for moderate 
or low Reynolds numbers, the effects of the numerical boundary layers can 
still be problematic 
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Non- conforming boundaries add an extra layer of 
difficulty. The search for means to better control these numerical artifacts is 
an active area of research. 

The numerical boundary layers in projection methods reflect in the known 



convergence results for them (e.g. |3ll l33|, |36|). Convergence is stated only 
in terms of integral norms, with the main difficulties near the boundary. 
There point-wise convergence (and even less convergence of the flow veloc- 
ity gradient) cannot be guaranteed — even if the solution is known to be 
smooth. Hence the accurate calculation of wall stresses with these methods 
is problematic. Guermond, Minev and Shen [l2| provide further details on 
convergence results, as well as an extensive review of projection methods and 
the improved pressure-correction schemes. 



Two other methods for solving the Navier-Stokes equations are the 



im- 



mersed boundary 24, 25, 29, 30, 37|, and the vortex-streamfunction (ii Ii. 26 
methods. These also decouple the calculation of the velocity and of the pres- 
sure. The immersed boundary method does so by introducing Dirac forces to 
replace the domain walls, which makes obtaining high order implementation 
of the boundary conditions difficult. The vortex-streamfunction formulation 
decouples the equations, but at the expense of introducing integral bound- 
ary conditions for the vorticity, which negatively impacts the efficiency. An 
interesting variation of the vortex-streamfunction approach, using only local 
boundary conditions, is presented in reference 
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Closely related to the immersed boundary methods are the penalty (al- 
ternatively: fictitious domain or domain embedding) methods — e.g. see [l|, 
[s, 19]. These methods, effectively, replace solid walls in the ffuid by a porous 
media with a small porosity < r/ <^ 1. In the limit r/ — )■ 0, this yields no slip 
and no ffow-through at the solid walls. Two important advantages of this 
approach are that complicated domains are easy to implement, and that the 
total fluid-solid force can be computed using a volume integral, rather than 
an integral over the boundary of the solid. Unfortunately, the parameter r) 
introduces boundary layers which make convergence slow and high accu- 
racy computations expensive, since rj cannot be selected independently of the 
numerical grid size. Spectral methods j^, 0] are very efficient and accurate, 
but they have geometry restrictions. 

Finally, we mention the algorithms based on a Pressure Poisson Equation 



fPPE) reformulation of the Navier-Stokes equations [ll|, |l6|, ll7|, |2l|, |32|, |3J, 



35| . which is the class of methods within which the work presented in this 
paper falls. In this approach the incompressibility constraint for the flow 
velocity is replaced by a Poissoiil equation for the pressure. This then allows 
an extra boundary condition — which must be selected so that, in fact, in- 
compressibility is maintained by the resulting system. This strategy was first 
proposed by Gresho and Shani , who pointed out that adding V ■ u = 
as a boundary condition yields a system of equations that is equivalent to 
the Navier-Stokes equations. Unfortunately, their particular PPE formula- 
tion incorporates no explicit boundary condition that can be used to recover 
the pressure from the velocity, by solving a Poisson problem — for a more 
detailed discussion of this, see remark H] in this paper. In FJ, Johnston and 



Liu propose a PPE system where this issue is resolved. In this paper — in 
equations fl2UH^ — we present another PPE system, also equivalent to the 
Navier-Stokes equations, which allows an explicit recovery of the pressure 
given the fiow velocity. A comparison with the one in [l7| can be found in 
remark [5] in this paper. 

PPE reformulations of the Navier-Stokes equations, such as the one in 
equations (plH^Ti) . or in reference [l^, have important advantages over the 
standard form of the equations. First, the pressure is not implicitly coupled 
to the velocity through the momentum equation and incompressibility. Hence 
it can be directly (and efficiently) recovered from the velocity field by solving 



^The choice of the Poisson equation for the pressure is not unique, e.g. see [If 
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a Poisson equation. This allows the velocity field to be marched in time using 
the momentum equation, with the pressure interpreted as some (complicated) 
function of the velocity. Second, no spurious boundary layers are generated 
for neither the velocity, nor the pressure. This because: 

— There are no ambiguities as to which boundary conditions to use for 
the pressure — hence errors induced by not-quite-correct boundary 
conditions do not occur. Note that these errors should not be confused 
with the truncation errors that any discretization of the equations will 
produce. Truncation errors are controlled by the order of the approxi- 
mation and, for smooth solutions, are uniformly small. 

— Incompressibility is enforced at all times. 

Hence pressure and velocities that are accurate everywhere can be obtained, 
in particular: near the boundaries. Finally, PPE formulations allow, at least 
in principle, for the systematic generation of higher order approximations. 

It follows that PPE strategies offer the promise of a resolution of the 
dilemma alluded to in the second paragraph of this introduction. They retain 
many of the advantages that have made projection methods popular, while 
not suffering from the presence of numerical boundary layers. On the negative 
side, the boundary conditions for PPE systems tend to be more complicated 
than the simple ones that the standard form of the equations has. 

This paper is organized as follows. In § [2] we discuss the Pressure Pois- 
son Equation (PPE) formulation of the Navier-Stokes equations. In § [3] we 
present a PPE formulation, using an alternative form of the boundary condi- 
tions, that allows a complete splitting of the momentum and pressure equa- 
tions. Namely: the pressure can be recovered from the flow velocity without 
boundary condition ambiguities. In §|l]we address a numerical question (sec- 
ular growth of the error under some conditions) and introduce a corrected 
formulation, involving a feedback parameter A, where this potential problem 
is corrected. This system, which is equivalent to the Navier-Stokes equa- 
tions, is displayed in equations ( 12UH2T]) . A discussion and comparison with 



the alternative PPE formulation by Johnston and Liu [17j is also included 
in this section. In § [5] we describe a second order solver for the new system 
introduced in § HI for irregular domains embedded within a cartesian stag- 
gered grid. In §[6], we implement and test our proposed scheme. This section 
includes convergence plots, indicating second order uniform convergence all 
the way to the boundary (L°° norm), of the pressure, the flow velocity, and 
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of the derivatives of the flow velocity. Finally, §[7] has the conclusions. The 



contents of the two appendices is as follows. In Appendix A we present an 



extended system of equations, valid for arbitrary flows, which for smooth 



solutions has the Navier-Stokes equations as an attractor. In Appendix B 
for completeness, we display formulas for the V ■ u = boundary condition 
for general conforming boundaries in curvilinear coordinates. 

2. The Pressure Poisson Equation. 

In this section we introduce the well-known pressure Poisson equation 
(PPE), and use it to construct a system of equations (and boundary con- 
ditions) equivalent to the const ant- density (hence incompressible) Navier- 
Stokes equations, with the velocity prescribed at the boundaries. Specifically, 
consider the incompressible Navier-Stokes equations in a connected domain 
n G R^, where D = 2 or D = 3, with a piece- wise smooth boundary 
Inside Q, the flow velocity field u(x, t) satisfies the equations 

ut + (u-V)u = /iAu-Vp + f, (1) 
V-u = 0, (2) 

where /i is the kinematic viscosityH p(x, t) is the pressure, f(x, t) are the 
body forces, V is the gradient, and A = is the Laplacian. Equation ([1]) 
follows from the conservation of momentum, while ([2]) is the incompressibility 
condition (conservation of mass). 

In addition, the following boundary conditions apply 



u 



X, t) for X G dn, (3) 
n-gdA = 0, (4) 



where 



Ian 

n is the outward unit normal on the boundary, and dA is the area (length 
in 2D) element on dQ. Equation (jl]) is the consistency condition for g, since 
an incompressible fluid must have zero net flux through the boundary. 



■^We work in non-dimensional variables, so that fi — 1/Re (where Re is the Reynolds 
number) and the fluid density is p = 1. 
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Finally, we assume that initial conditions are given 

u(x, 0) = uo for X e 11, (5) 

V ■ uo = for X G 1), (6) 

uo(x) = g(x, 0) for X e dn. (7) 



Remark 1. Of particular interest is the case of fixed impermeable walls, 
where no flux u ■ n = and no slip u x n = apply at dQ. This corresponds 
to g = in ([3]). Note that the no-slip condition is equivalent to u ■ t = for 
all unit tangent vectors t to the boundary. Jit 

Remark 2. In this paper we will assume that the domain Q is fixed. Situa- 
tions where the boundary of the domain, dfl, can move — either by externally 
prescribed factors, or from interactions with the fluid, are of great physical 
interest. However, to keep the presentation of the new method as simple as 
possible, we postpone consideration of these cases for further work. Jft 

Next we introduce the Pressure Poisson Equation (PPE). To obtain the pres- 
sure equation, take the divergence of the momentum equation ([T]), and apply 
equation (|2]) to eliminate the viscous term and the term with a time deriva- 
tive. This yields the following Poisson equation for the pressure 

Ap = V ■ (f - (u ■ V) u) . (8) 

Two crucial questions are now (for simplicity, assume solutions that are 
smooth all the way up to the boundary) 

[2^ Can this equation be used to replace the incompressibility condition ^? 
Since — given ^ — the divergence of ([T]) yields the heat equation 

0t = /^A</., (9) 

for = V ■ u and x G fi, it would seem that the answer to this 
question is yes — provided that the initial conditions are incompressible 
{i.e.: = for t = 0). However, this works only if we can guarantee 
that, at all times, 

= 0, (10) 

for X G dVt. 
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Given the flow velocity u, can ^ be used to obtain the pressure p? 
Again, at first sight, the answer to this question appears to be yes. 
After all, ([8]) is a Poisson equation for p, which should determine it 
uniquely — given appropriate boundary conditions. The problem is: 
what boundary conditions? Evaluation of ([1]) at the boundary, with 
use of ([3]), shows that the fiow velocity determines the whole gradient 
of the pressure at the boundary, which is too much for ([H]). Further, if 
only a portion of these boundary conditions are enforced when solving 
dH]) — say, the normal component of ([1]) at the boundary, then how can 
one be sure that the whole of ([1]) applies at the boundary? 

Remark 3. From an algorithmic point of view, an affirmative answer to the 
questions above would very useful, for then one could think of the pressure 
as some (global) function of the flow velocity, in which case ([1]) becomes 
an evolution equation for u, which could then be solved with a numerical 
"marching" method. ^ 

The issue in item|2^ can be resolved easily, and we do so next. We postpone 
dealing with the issue in item |2}3 till the next section, § [31 Since the addition 
of an equation for the pressure allows the introduction of one extra bound- 
ary condition, we propose to replace the system in ([T]), ([2]), and by the 
following Pressure Poisson Equation (PPE) formulation: 

Ut + (u-V)u = /iAu-Vp + f, (11) 
Ap = V ■ (f - (u ■ V) u) , (12) 

for X G f2, with the boundary conditions 

u = g(x, t), (13) 
V-u = 0, (14) 

for x G dQ — where, of course, the restriction in (jl]) still applies. The extra 
boundary condition is precisely what is needed to ensure that the pressure 
enforces incompressibility throughout the flow (see item [2^) 

Remark 4. It can be seen that for smooth enough solutions (u, p), the pair 
of equations f lTT| - [T2|) . accompanied by the boundary conditions f|T3] - [T^ . are 
equivalent to the incompressible Navier-Stokes equations in ([1]), ([2]), and 
(|3]). Of course, this result is not new. This reformulation of the Navier- 



Stokes equations was flrst presented by Gresho and Sani [ll[ — it can also 
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be found in reference 
Welch 
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32| . However, it should be pointed out that Harlow and 
in their pioneering work, had already noticed that the boundary 
condition V ■ u = was needed to guarantee, within the context of their 
MAC scheme, that V ■ u = everywhere. 

This formulation does not provide any boundary conditions for the pres- 
sure, which means that one ends up with a "global" constraint on the so- 
lutions to the Poisson equation f[T^ . Hence it does not yield a satisfactory 
answer to the issue in item [2]d, since recovering the pressure from the flow 
velocity is a hard problem with this approach. 

Direct implementations of (ITTHT^ have only been proposed for simple ge- 
ometries — in particular: grid-conforming boundaries. In 2l| a spectral al- 
gorithm for plane channel flows is presented. We have already mentioned 14| . 
where they use a staggered grid on a rectangular domain, and the condition 
V • u = is used (at the discrete level) to close the linear system for the 



pressure. In [16|, by manipulating the discretization of the boundary con- 
ditions and of the momentum equation ffTTl) . they manage to obtain "local" 
approximate Neumann conditions for the pressure — both on rectangular, 
as well as circular, domains where u = on the boundary. Unfortunately, 
the approaches in [l6[ seem to be very tied up to the details of particular 
discretizations, and require a conforming boundary. X 



3. Theoretical Reformulation 

We still need to deal with the issue raised in item In particular, in 
order to implement the ideas in remark [3l we need to split the boundary 
conditions for (|TT1 - [T2l) . in such a way that: (i) there is a specific part of the 
boundary conditions that is used with f lTT]) to advance the velocity field in 
time, given the pressure, (ii) The remainder of the boundary conditions is 
used with f[T^ to solve for the pressure — at each fixed time, given the flow 
velocity field u. This is the objective of this section. 

The conventional approach in projection or fractional step methods is to 
associate (fT3|) with equation ( ITT]) . This is reasonable for evolving the heat- 
like equation (fTT]) . Unfortunately, it has the drawback of only implicitly 
defining the boundary conditions for the Poisson equation ( IT^ . Specifically, 
the correct pressure boundary conditions are those that guarantee V ■ u = 
for X G dQ, and such boundary conditions cannot easily be known a priori 
as a function of u. Hence one is left with a situation where the appropriate 
boundary conditions for the pressure are not known. This leads to errors 
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in the pressure, and in the incompressibihty condition, which are difficult to 
control. In particular, errors are often most pronounced near the boundary 
where no local error estimates can be produced, even for smooth solutions. 
By the latter we mean that the resulting schemes cannot be shown to be 
consistent, all the way up to the boundary, in the classical sense of finite 



differences introduced by Lax [23 



In this paper we take a different approach, which has a similar spirit to the 



one used by Johnston and Liu [17| — see remark [51 Rather than associate 
all the D components of (fT3|) with equation ( [TT|) . we enforce the D — 1 
tangential components only, and complete the set of boundary conditions for 
(11 ip with (I14p . Hence, when evolving equation (ITT]) we do not specify the 
normal velocity on the boundary, but — through the divergence condition 
f lT^ — specify the normal derivative of the normal velocity. Finally, the 
(as yet unused) boundary condition on the normal velocity, n ■ (u — g) = 
for X G dVl, is employed to obtain an explicit boundary condition for 
equation f lT2|) . We do this by requiring that the pressure boundary condition 
be equivalent to (n ■ (u — g))^ = for x G dVL — which then guarantees that 
the normal component of f[T^ holds, as long as the initial conditions satisfy 
it. This objective is easily achieved: dotting equation (fTTj) through with n, 
and evaluating at the boundary yields the desired condition. The equations, 
with their appropriate boundary conditions, are thus: 

Ut — /i Au = — Vp — (u ■ V) u + f for X G 
n X (u - g) = for X G )- (15) 

V ■ u = for X G 911, 



and 



Ap = -V-((u-V)u) + V-f for xG n, 

n ■ Vp = n ■ (f — gt + yU Au — (u ■ V) u) for x G d^l. 



(16) 



Again, for smooth (up to the boundary) enough solutions (u, p) of the equa- 
tions: the incompressible Navier-Stokes equations ([IH2]), with boundary con- 
ditions as in ([3]), are equivalent to the system of equations and boundary 
conditions in f lTSHTB]) . 

For the sake of completeness, we display now the calculation showing 
that the boundary condition splitting in ( [T5] - [l6|) recovers the normal velocity 
boundary condition n ■ u = n ■ g. To start, dot the first equation in ( [T5|) with 
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the normal n, and evaluate at the boundary. This yields 

n ■ = n ■ (/i Au - Vj9 - (u ■ V) u + f ) for x G dQ. (17) 

Next, eliminate n ■ Vp from this last equation — by using the boundary 
condition for the pressure in f|T6|) . to obtain 

(n ■ (u - g))^ = for X G dQ. (18) 

This is a trivial ODE for the normal component of the velocity at each point 
in the boundary. Thus, provided that n ■ (u — g) = initially, it holds for all 
time. 

An important final point to check is the solvability condition for the 
pressure problem. Given the flow velocity u at time t, equation (fT6|) is a 
Poisson problem with Neumann boundary conditions for the pressure. This 
problem has a solution (unique up to an additive constant) if and only if the 
"flux equals source" criteria 

/ n- (f-gi + /iAu- (u- V)u) dA= /" V- (f- (u- V)u) dV^ (19) 
Jan Jn 

applies. This is satisfied because: 
[li /^^ n ■ (f - (u ■ V) u)dA = V ■ (f - (u ■ V) u) d^, 

m /,^n-AudA = ^A(V-u)d\/ = 0, 
Et ^^n-g,dA = | /,^n-gdA = 0, 
where we have used Gauss' theorem, incompressibility, and (j4]). 

4. Modification for Stability 

For the system of equations in f|T5] - [T6l) . it is important to notice that 

UK The tangential boundary condition on the flow velocity, n x (u — g) =0 
for X G dQ, is enforced explicitly. 

IDd The incompressibility condition, V ■ u = for x G fi, is enforced 
"exponentially". By this we mean that any errors in satisfying the 
incompressibility condition are rapidly damped, because = V ■ u 
satisfies flO HTOj) . Thus this condition is enforced in a robust way, and 
we do not expect it to cause any trouble for "reasonable" numerical 
discretizations of the equations. 
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HI: By contrast, the normal boundary condition on the flow velocity, n • 
(u — g) = for X G dQ, is enforced in a rather weak fashion. By 
this we mean that errors in satisfying this condition are not damped at 
all by equation f|T8|) . Thus, this condition lacks the inherent stability 
provided by the heat equation. In practice, numerical errors add (effec- 
tively) noise to equation flTK]) . resulting in a drift of the normal velocity 
component. This can have de-stabilizing effects on the behavior of a 
numerical scheme. Hence it is a problem that must be corrected. 

In this section we alter the PPE equations f lTSHTBl) to address the problem 
pointed out in item Ht. We do this by adding an appropriate "stabilizing" 
term. Our goal here is to develop a pair of differential equations — fully 
equivalent to (fT5l - [T6|) — which are suitable for numerical implementation. In 
order to resolve the issue in itemH}:;, we add a feedback term to the equations, 
by altering the pressure boundary condition. Specifically, we modify the 
equations from (ITBHTB]) to 

ut — fiAu = — Vp— (u-V)u + f for xG ^l, 

n X (u - g) = for X G ). (20) 

V ■ u = ior xedn, 

and 

Ap = -V ■ ((u ■ V) u) + V ■ f for X G fi, 

n ■ Vp = n ■ (f - gt + /i Au - (u • V) u) } (21) 

+ A n ■ (u — g) for X G dQ, 

where A > is a numerical parameter — see § 14.11 This system is still 
equivalent to the incompressible Navier-Stokes equations and boundary con- 
ditions in (IIH3]), for smooth (up to the boundary) enough solutions (u, p). 
The heat equation fl9]- [T0l) for = V ■ u still applies, while the equation for 
the evolution of the normal velocity at the boundary changes from f|T8l) to: 

(n- (u-g))^ = -An- (u-g) for x G 91]. (22) 

Thus, if n ■ (u — g) =0 initially, it remains so for all times. In addition, this 
last equation shows that this new system resolves the issue pointed out in 
itemHt. 
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Finally, we check what happens to the solvability condition for the pres- 
sure problem, given the system change above. Clearly, all we need to do is 
to modify equation f|T9|) by adding — to its left hand side, the term 

\ [ n- (u-g) dA = \ [ V -udV - \ [ n-gdA = 0, (23) 

JdQ Jq Jdn 

where we have used incompressibility, and equation It follows that solv- 
ability remains valid. 

Remark 5. The reformulation of the Navier-Stokes equations in (I20H2T]) is 



similar to the one used by Johnston and Liu in [17[. In their paper the 



authors propose methods of solution to the Navier-Stokes equations based 



on the equivalent system (for simplicity, we set g = 0, as done in 17|) where 

(a) For X G Q, the same equations as in fl20] - [2T|) apply. 

(b) For X G dQ, u = is used for the momentum equation. 

(c) For X G dQ, n ■ Vp = n ■ (— /i V x V x u + f ) is used for the Poisson 
equation. This is equivalent to n ■ Vp = n ■ (/i A u — yU V(V ■ u) + f ). 

For this system the incompressibility condition = V ■ u = follows because 
these equations yield 

(d) 0f = ;U A for X G fi, with n ■ V = for x G d^l. 

Thus, if vanishes initially, it will vanish for all times. 

(e) The main advantage of this reformulation over the one in fl20] - [2Tl) is 
that (in general) the boundary condition V • u = in fl2U]) couples all 
the components of the flow velocity field — see § 15.31 Thus an implicit 
treatment of the viscous terms in (120] - [2T!) would be more expensive 
(and complicated) than for the system in items a-c above. However, it 
is not clear to us at this moment how much of a problem this is. The 
reason is that the coupling is "weak", by which we mean that: in the 
Ng X Ng discretization matrix for the Laplacian — where Nq is the 
number of points in the numerical grid, the coupling induced by the 

1 /2 1 /3 

boundary condition affects only 0{NJ ) entries in 2-D, and 0{NJ ) 
entries in 3-D — at least with the type of discretization that we use in 
§ 15.31 Hence, it may be possible to design algorithms where the extra 
cost is moderate, and not a show-stopper. 
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(f) The main disadvantage of this reformulation over the one in fl20] - [2Tl) is 
related to the discussion in items HJd and Hl^ earlier in this section — 
which lead us to add the feedback control parameter A in fEUH^T]) . In 
other words, the incompressibility condition is enforced in a relatively 
weak fashion: the heat equation evolution in item d forces errors in 
to converge (exponentially) to their mean, but if the mean of the 
errors is not zero, will drift away from zero. As shown by our nu- 
merical experiments in § |6l drift can be a problem (if not controlled 
by a parameter such as A) when doing calculations in domains with 
non-conforming boundaries. 

For simple geometries with conforming boundaries this is not a prob- 



lem, as shown by the numerical experiments in [17[. Unfortunately, 
they did not test their method with non-conforming boundaries, hence 
we do not actually know how much of a problem this can be. At any 
rate, if it is a problem, it is possible to modify the method by adding 
a parameter that forces exponential decay of to zero. For example: 
modify the Poisson equation for the pressure to 

Ap = -V ■ ((u ■ V) u) + V ■ f + A V ■ u. 

This then changes the system in item d to 

0t = /iA0 — A0 for X G fi, with n ■ V = for x G dQ. 



Finally, a minor disadvantage of the system in items a-c, relative to that in 
( 120] - [2T1) . is that: ( 120] - [2T|) reduces to the steady state Navier-Stokes equations 
for steady states, while (a-c) does not — since it allows V-u = constant ^ 0. 
This is, of course, a manifestation of the same issue raised in item f. X 

Remark 6. It would be nice to be able to use A = A(x), so as to optimize 
the implementation of the condition n ■ u = n • g for different points along 
dQ. However, this is not a trivial extension, since it destroys the solvabil- 
ity condition for the pressure. Thus, other (compensating) corrections are 
needed as well. We postpone the study of this issue for future work. X 

Remark 7. As show earlier in f lT9|) and (1231) . the validity of the solvability 
condition, for the problem in (^T^, relies on the velocity field satisfying the 
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incompressibility constraint V ■ u = 0. On the other hand, in the course of 
a numerical calculation, the discretization errors result in a small, but non- 
zero, divergence — thus solvability fails. However, the errors in solvability 
are small. Hence, a least squares solution of the discretized linear equation 
for the pressure provides an approximation within the order of the method 
— which is as good as can be expected. These considerations motivate the 
following theoretical question: can the equations in l[W^WJ\) be modified, so 



they make sense even for V ■ u 7^ ? In Appendix A we show that this is 



possible. X 

4.I. Selection of the parameter A. 

For numerical purposes, here we address the issue of how large A should 
be, by using a simple model for the flow's normal velocity drift. Notice that 
no precision is needed for this calculation, just order of magnitude. In actual 
practice, one can monitor how well the normal velocity satisfies the boundary 
condition, and increase A if needed. In principle one should be wary of using 
large values for A, since this will yield stiff behavior in time. However: (i) the 
(viscous) momentum equation in ( 120|) is already very stiff. Hence, A would 
have to be fairly large before it increases overall stiffness, (ii) The calculation 
below shows that A does not need to be very large. 

It seems reasonable to assume that one can model how the numerical 
errors affect the ODE fl22l) for £ = n ■ (u — g), by perturbing the coefficients 
of the equation, and adding a forcing term to it. Hence we modify equation 
(|22|) as follows 

£^ = -\cp£ + e-f for X G dn, (24) 

where e -C 1 characterizes the size of the errors (determined by the order of 
the numerical method), while Cp = Cp(x, t) = l+0{e) and 7 = 7(x, t) = 0(1) 
are functions encoding the numerical errors. What exactly they are depends 
on the details of the numerical discretization, but for this calculation we do 
not need to know these details. All we need is tha10 

< Cm < Cp and I7I < T, (25) 

where Cm and T are some positive constants, with Cm ~ 1, and T = 0(1) 
— but not necessarily close to one. 



^For smooth enough solutions, where the truneation errors are controhcd by some 
derivative of the solution. 
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The solution to is given by 



£ = £oe-^^' +e / 7(x, s) e'^^^ ds, (26) 



J 



where £o is the initial value, Ji = /i(x, t) = Cp(x, s) ds, and I2 = /i(x, t) — 
Ji(x, s). The crucial term is J, since the first term decreases in size, and 
starts at the initial value. However 

\J\<V l\-^^^'^'~'^ds<^—. (27) 

Within the framework of a numerical scheme, the normal boundary ve- 
locity may deviate from the prescribed velocity by some acceptable error 8. 
Thus we require e J = 0(5) or less, which — given ( 1271) — will be satisfied if 

6 r 6 r 

A ~ ~ -T' larger. (28) 

For the second order numerical scheme in § [5l it is reasonable to expect that 
e = (Ax)2, and to require that 5 = {^xf . Then (EH]) reduce^H to A > T. Of 
course, we do not know (a priori) what F is; this is something that we need 
to find by numerical experimentation — see the first paragraph in this § 14.11 
For the numerical calculations reported in § El we found that values in the 
range 10 < A < 100 gave good results. 

5. Numerical Scheme 

In this section we outline an efficient numerical scheme for solving the 
coupled differential equations ( |20l - [2T]) on a two-dimensional irregular domain. 
As should be clear from the prior sections, the main issue we aim to address 
in this paper, is that of how to effectively implement the incompressibility 
condition and the boundary conditions for the pressure, avoiding the difficul- 
ties that projection and fractional step methods have. These are problems 
that are not related to the nonlinearities in the Navier-Stokes equations, and 
occur even for the linearized equations. Hence, in the spirit of focusing on 



^Note that this estimate for A is independent of the grid mesh-size Aa;. 
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the key issues only — also see remark [SI the calculations presented in § 
are for the linear equations onlyl^l Furthermore, while in the description of 
the scheme in this section we carry through the nonlinear terms, we do not 
describe their implementation. This should not be interpreted as implying 
that we believe that this is a trivial matter. Quite the opposite. However, 
there is an extensive literature on this topic 0] and we have nothing to add 
to it in this paper. 

To achieve an efficient scheme, we decouple the pressure and velocity 
fields, and explicitly treat each term in the time evolution of ( 120] - [2T|) . Specif- 
ically, since both the right hand side and boundary conditions of equation 
(PTj) depend solely on u, we may view the pressure as a computable functional 
of the velocity, p = p[u]. The computation of p[u] requires the solution of 
a Poisson equation with Neumann boundary conditions. With this in mind, 
the momentum equation then has the form = F[u], where F[u] has a 
complicated, yet numerically computable form: 

F[u] = Au - Vp[u] - (u ■ V) u + f for X en. (29) 

We now use an explicit forward Euler scheme to discretize the time evolution 
for fl20] - [2T|) . paired with an appropriate discretization in space described later 
in this section. This yields the scheme 

^ (u"+^ - u") = /i Au" - Vj9" - (u" • V) u" + f " for ^ e Q, (30) 

with boundary conditions n x u = n x g and V ■ u = for x G dQ, where 
the pressure is given by 

Ap" = -V ■ ((u'^ ■ V) u") + V ■ r for X G f], (31) 

with the boundary condition 

n • Vp" = n ■ (f" - + /i Au" - (u" ■ V) u") + An ■ (u" - g") 

for X G dQ. Here, starting with the initial data u", a superscript n is used 
to denote a variable at time t = n At, where < At ^ 1 is the time step. 

^Though we have carried some calculations with the nonlinear terms, these are not 
presented here. 
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Remark 8. Our purpose here is to illustrate the new approach with a sim- 
ple scheme that does not obscure the ideas in the method with technical 
complications. Hence, the scheme here is first order in time (explicit) and 
second order in space, with the stability restriction At oc (Ax)^. However, 
unlike projection methods and other approaches commonly used to solve the 
Navier-Stokes equations — see § [H this new formulation does not seem to 
have any inherent order limitations. Unfortunately, the Navier-Stokes equa- 
tions are stiff and nonlinear, which means that the fact that higher order 
extensions are possible does not mean that they are trivial. X 

5.1. Space grid and discretization. 

To discretize the equations in space, we use finite differences over a carte- 
sian, square (Ax = Ay), staggered grid. The pressure values are stored at 
the nodes of the grid, while the horizontal and vertical components of the 
velocity are stored at the mid-points of the edges connecting the grid nodes 
(horizontal component on the horizontal edges and vertical component on 
the vertical edges). 

When handling an arbitrary curved boundary, we cannot conform the 
boundary to the grid, but rather we immerse it within the regular mesh — 
see figure [1] Then, to numerically describe the domain boundary, we identify 
a set Ch of Ne points in dfl, say x^j = {xb, Vhjj for 1 < j < — see itemEt 
below. These Nf, points are located at 0{Ax) distances apart, so that the 
resolution of the boundary is comparable with that of the numerical grid. 

For the heat equation Tt = fi AT, the stability restriction for the standard 
scheme using a 5 point centered differences approximation for the Laplacian, 
and forward Euler in time, is 

At<C^-^, (32) 

where C = and D = 1, 2, ... is the space dimension. Since the method 
described here uses exactly the same approach to advance the velocity flow 
field u, we expect the same restriction (with, perhaps, a different constant C) 
to apply. For the 2D numerical calculations presented in § El we found that 
the algorithm was stable with C < 0.2, while C > 0.3 generally produced un- 
stable behavior. Below we separately address the numerical implementation 
of equations (130|) and (!3T|) . 
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Figure 1: This plot shows the staggered grid, and the boundary. The numer- 
ical pressure values correspond to the graph nodes, while the velocities (ar- 
rows) correspond to the edge midpoints. Here the circles (o) and squares (□) 
denote ghost pressure points, and boundary velocities, respectively. These 
are used to implement the boundary conditions in the Poisson and momen- 
tum equations, respectively. The diamonds (0) denote the points {xb, yb)i, 
used to represent the boundary dQ. 



5.2. Poisson Equation 

To solve the pressure Poisson equation (I3T!) with Neumann boundary 



to 



conditions, we use the ghost point idea with a method by Greenspan [10 
implement the boundary conditions. As discussed above, we embed the do- 
main Q within a cartesian square grid, and classify the computational points 
in the grid into inner and ghost points (see next paragraph). Then: (i) for the 
inner points we discretize the Laplace operator using the standard 5 point 
centered differences stencil, (ii) For the ghost points we obtain equations 
from an appropriate discretization of the Neumann boundary condition. In 
particular, if there are Na inner points, and Nb ghost points, then the dis- 
cretized pressure is represented by a vector p G where N = Na + Nb. 

The computational domain, Cp, for the pressure is comprised by the inner 
points and the ghost points, defined as follows: 

The inner points are the points in the cartesian grid located inside fl. 
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\5jp The ghost points are the points in the cartesian grid that are either 
outside Q or on dQ, and which are needed to complete the 5 point 
stencil for some inner point. 

Construct the set Cb used to track the boundary dfl as follows: For 
every pressure ghost point in the grid: x^^, 1 < j < Nf,, select the 
closest point in the boundary x^j E dQ. We will use these points to 
produce equations that approximate the Neumann boundary conditions 
(one equation per point). 

Remark 9. With the set Cb as above, Nf, = Nb- This is not the only pos- 
sibility. In the end, the discretized equations for the pressure have to be 
solved in the least squares sense (see remark [10]). Thus, we could have 
elected to over-determine the implementation of the boundary condition by 
selecting Ne > Nb points in dQ. One reason for doing this being to obtain 
a "smoother" implementation, avoiding the irregularities that occur as the 
boundary placement (relative to the cartesian grid) changes along dQ. In 
our calculations we found that this is not a must for the solution of the Pois- 
son equation. On the other hand, for the momentum equation this strategy 
is needed to avoid numerical instabilities — see § 15.31 Jit 

It follows that, on the computational domain Cp, the Poisson equation 
(including the boundary conditions) is discretized by Na equations inside fl, 
and Nb equations derived from the boundary conditions. Specifically 

Oi We use the standard 5 point centered differences stencil for the Lapla- 
cian, to obtain one equation for each of the Na inner points. Similarly, 
we use a second order approximation for the forcing terms on the right 
hand side of the equation. 

We construct one boundary condition equation for each pressure ghost 
point X™, 1 < j < Nb. This is done by using 6 nearby points from the 
computational domain to obtain a second order approximation to the 
normal derivative at the corresponding point x^j E dQ — see item [St 
above. We use a similar approach to approximate the terms on the 
right hand side of the Neumann boundary conditions. 

Hence, each boundary condition equation involves the ghost point, 
points inside the domain, and (possibly) other nearby ghost points. 
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Thus the Poisson equation can be written in a matrix form with the following 
structure 

Here B G R^''^^ and L G R^"^^ are the discrete matrix representations for 
the derivatives (n ■ V) in the Neumann boundary condition,, and Laplacian 
(A), respectively. The terms a G R^" and b G H^'' are the discrete represen- 
tations of the source terms and applied boundary conditions for the Poisson 
equation: 

a ^ V- (f- (u- V)u), (34) 
b ^ n- (f-gt + /iAu- (u- V)u) + 

An-(u-g). (35) 

Thus the pair of equations 

Lp = a and Bp = b (36) 
are the discrete analog of equation fl2Tl) . 

Remark 10. As pointed out earlier in remark [71 the discretization errors 
cause solvability for equation fl33p to fail. These solvability errors are "small" 
(second order). Hence, as it is commonly the case with numerical solutions 
of the Poisson equation with Neumann boundary conditions, we solve equa- 
tion fl33|) in the least squares sense. This provides an approximation to the 
pressure which is within the discretization error. 4> 

5.3. Momentum Equation 

The main numerical difficulty with implementing f l3Up is produced by 
the boundary conditions. Specifically: on a cartesian staggered grid the 
boundary conditions V ■ u = and n x (u — g) = couple the "horizontal" , 
u, and "vertical", v, components of the fiow velocity — with the exception 
of the special case of a boundary aligned with the grid, where there is no 
coupling. Hence, in general, the implementation of the boundary conditions 
requires the solution (at every time step) of a linear system of equations that 
couples all the boundary velocities (these are defined below). 

The computational domain for the velocities, Cu, is defined in terms of 
the edges in the cartesian grid with which the velocities are associated (see 
figured]). To define Cu, it is convenient to first introduce the extended set of 
pressure nodes, Sp, from which is easily constructed: 



21 



Ef A pressure node in the cartesian grid belongs to £p if and only if it 
either belongs to Cp, or if it is connected (by an edge in the grid) to a 
ghost pressure node that lies on dVL. 

A velocity is in Cu if and only if: (a) Its corresponding edge connects 
two points in £p. (b) At least one of the two points is in Vt (or dVt). 

Remark 11. Notice that £p is exactly what Cp becomes if dVl is modified 
by an "infinitesimal" perturbation that turns all the ghost pressure points 
on dVl into inner pressure points. 

The computational domain is, in turn, sub-divided into inner and bound- 
ary velocity edges 

[5)1 An edge in Cu is an inner velocity if and only if Cu includes the four 
other edges needed to compute the Laplacian (either Am or Aw) at the 
edge mid-point, using the 5 point centered differences approximation. 

[5} An edge in is a boundary velocity if and only if it is not an inner 
velocity. Let M he the number of boundary velocities. 

The solution of the momentum equation (!20l) is thus performed as follows: 

[3] At the start of each time step the right hand side in (1501) is approxi- 
mated (to second order, using centered differences) at the inner velocity 
locations, which can then be updated to their values at the next time. 

(Sk Next, using the boundary conditions, the values of the boundary ve- 
locities at the next time step are constructed from the inner velocities. 
This "extension" process is explained below. 

Let y G R*^ be the vector of boundary velocities. Then y is determined by 
two sets of equations, corresponding to the discretization on dVL of V ■ u = 
and n X (u — g) = 0. In our approach the divergence free criteria is enforced 
"point-wise", while the condition on the tangential velocity is imposed in a 
least squares sense. 

Implementation of the divergence free V ■ u = boundary condition. At 
first sight, this boundary condition appears to be the hardest to implement, 
since it is a Neumann condition (essentially) prescribing the value of the 
normal derivative of the normal component of the velocity, in terms of the 
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tangential derivative of the tangential velocity. However, because (for the 
exact solution) V ■ u = everywhere, an implementation of this condition 
which is second order consistent (in the classical sense of finite differences 



introduced by Lax [23|) is easy to obtain, as follows: 

First: identify the Md pressure nodes in Cp, which have at least one 
boundary velocity as an adjacent edge. These pressure nodes lie either 
on (9f2, or inside fi, and are all within a distance 0{Ax) of dVL — definitely no 
further away than a/2 Ax. Second: for each of these Md points, use centered 
differences to approximate the flow divergence at the point, and set V-u = 0. 
This provides Md equations that couple the (unknown) boundary velocities 
to the (known) inner velocities. In matrix form this can be written as 

Dy = s (37) 

where D is the portion of the discrete divergence operator acting on the 
unknown boundary velocities, while s is the associated flux derived from the 
known inner velocities. D is a rectangular, very sparse, matrix whose entries 
are and ±1 — note that Ax can be eliminated from these equations. 

Remark 12. Notice that V • u = for the exact solution. Hence, setting 
V ■ u = at the nodes near the boundary (as done above) involves no 
approximation. The second order nature of (!37|) is caused by the error in 
computing V-u — there is no "extrapolation" error. 



Implementation of the tangential velocity nx (u — g) =0 boundary condition. 
It is easy to see that Md ~ M/2, since most of the Md pressure nodes in the 
selected set will connect with two boundary boundary velocities. Thus f l37|) 
above provides approximately one half the number of equations needed to 
recover y from the inner velocities. It would thus seem natural to seek for 
M — Md additional equations using the other boundary condition. Namely, 
find M—Md points on dVL, and at each one of them write an approximation to 
n X (u — g) = using nearby points in Cu. Unfortunately, this does not work. 
It is very hard to do the needed approximations in a fashion that is robust 
relative to the way VL is embedded in the rectangular grid. Our attempts 
at this simple approach almost always lead to situations where somewhere 
along dVt an instability was triggered. 

To avoid the problem stated in the previous paragraph, we over-determine 
the implementation of the tangential velocity boundary condition, and solve 
the resulting system in the least squares sense. The boundary condition is 
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replaced by the minimization of a (discrete version) of a functional of the 
form 

/ \n X {u - g)f w dA, (38) 
Jan 

where w is some (strictly positive) weight function. This approach yields a 
robust, numerically stable, approximation — fairly insensitive to the partic- 
ular details of how Q is embedded within the cartesian grid. 
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Figure 2: This plot illustrates the implementation of the momentum equation 
boundary conditions. The circles (o) indicate the points at which the velocity 
divergence is set to zero. The six arrows indicate the horizontal velocity 
components in the patch, V'^^, used to extrapolate u to the three boundary 
points indicated by the diamonds {()). The squares (□) denote the boundary 
velocities — which are part of the patch VJ'^. 

In fact, we do not implement the minimization of a functional as in (138 p . 
but a simpler process which is (essentially) equivalent. To be specific, we 
start with the set Cb of points in dQ — see item [St. For each x^j G Cb, where 
1 < J < ^e, we identify a local horizontal, P", and vertical, VJ, velocity 
"patch" . These patches — see figure [2] — have the following properties 

[5j Each patch contains both inner and boundary velocities. 
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Ehi Each patch contains 6 velocities, in an appropriate structure, so that 
these velocities can be used to extrapolate (with second order accuracy) 
values of the corresponding velocity {u or v) to nearby points along the 
boundary. For example: the 5 point stencil used to approximate the 
Laplacian, plus one of the four next closest points to the center point. 

Eh The union of all the patches contains all the boundary velocities. 

These patches are used to (linearly) interpolate/extrapolate the velocities to 
nearby points along the boundary. Each patch is used to extrapolate the 
velocity to the three "closest" boundary points to the patch. In this fashion, 
for every x^j, 1 < j < A^^e, three different approximations to the tangential 
velocity at the boundary follow. These involve different (linear) combinations 
of velocities at the nearby edges in C^, including both boundary and inner 
velocities. In this fashion, a set of 3 Ng additional linear equations — beyond 
those in equation fl37|) — for the boundary velocities follow. 

Remark 13. The idea in the process above is to "link" the patches used to 
extrapolate the velocity to the boundary. This is so that no "gaps" occur in 
the implementation of the tangential boundary conditions, as the positioning 
of (9f2 relative to the grid changes. We elected to use minimal coupling 
(each patch linked to its two neighbors). In principle one could increase the 
amount of over-determination of the tangential boundary condition. This 
would require using a set of points in dfl that is "denser" than Cb, but in our 
calculations we found that this was not needed. 

Finally, an interesting question is: why is a process similar to this one not 
needed for the implementation of the boundary conditions for the Poisson 
equation in § 15.21 .'' The obvious answer is that the Poisson equation itself 
couples everything, so that no extra coupling needs to be added. X 

Putting everything together, an overdetermined system of equations for 
the boundary velocity vector y follows, which can be written in the form 



Here E G j^^AfexM^ ^^^q second equation is to be solved in the least 
squares sense, subject to the constraint imposed by the first equation. One 
way to do this is as follows: write 



Ey 



= s, 
= t. 



(39) 
(40) 



y = yp + 



(41) 
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where P is a matrix whose columns are a basis for the kernel of D — hence 
D P = 0, y-p is a particular solution of (15^ . and c is a constant vector 
parameterizing the space of (numerical) divergence free boundary velocities. 
Substituting the ansatz fHTl) into fHOl) yields 

{EP)c = t-Eyp, (42) 

which is a constraint-free least squares problem for the vector c. Thus we 
can write 

c={EPy (t-Ey,), (43) 
which, together with 0411) . gives the solution y. 

Remark 14. In fH3|) . {E P)^ is the pseudo-inverse, defined by the singular 
value decomposition of the matrix F = EP. In fact, since F'^ F should be 
invertible (see below), F+ = {F'^ F)'^ F'^ . 

One of the objectives of the above construction is to ensure that the 
boundary velocities are completely determined by the boundary conditions 
(and the inner velocities). If the problem in ( HOj) for y, with the constraint 
in (139 p . were to have more than one solution that minimizes the norm of 
t — Ey, then this would be a sure sign that the boundary condition imple- 
mentation is flawed, and there is missing information (i.e.: more equations 
are needed, see remark [T3|) . 

Of course, the requirement in the prior paragraph is equivalent to the 
statement that F is invertible. Finally, notice that in a domain with a 
fixed boundary, one may preprocess the matrices E, D and P. Jit 

Remark 15. It should be clear that the scheme developed in this section, 
is second order consistent0 (all the way up to the boundary) in the classical 



sense of finite differences introduced by Lax [23 



5.4- Comparison with the Projection Method 

The scheme proposed here — which involves the implementation of the 
equations in f l30] - [3T|) . appears to be quite similar to fractional step methods, 
such as Chorin's |6| original projection method. Specifically: advancing one 
time step requires both the evolution of a diffusion equation for the flow ve- 
locity, and the inversion of a Poisson equation for the pressure — which is the 



''For the system of PDE in 
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same as in projection methods. But there are important differences, mainly 
related to the implementation of the boundary conditions, and of the in- 
compressibility condition. Below we highlight the similarities and differences 
between the pressure Poisson approach proposed here and, for simplicity, the 
projection method in its original formulation |6i]. A (fairly recent) thorough 
review of projection methods, exploring the improvements to the approach 
since joj, as well as their drawbacks, can be found in 

First, the starting point for the method here is the discretization of a re- 
formulation of the equations: (120] - [2T]) . not the "original" Navier-Stokes equa- 
tions ([TH2])- In this equivalent set: (i) there is a natural way to recover the 
pressure from the flow field at any given time, (ii) The time evolution auto- 
matically enforces incompressibility. As a consequence, there is (in-principle) 
no limitation to the order in time to which the reformulated equations can 
be numerically discretized. In contrast, the pro ject ion method is equivalent 



to an approximate L U matrix factorization [12|, |28|, |37[ of the discrete differ- 



ential operators coupling u and p. This approximate factorization yields a 
splitting error in time, which is very hard to circumvent in order to achieve 
higher accuracy. 

Second, the pressure Poisson formulation used here ensures both that, at 
every point in the time integration: (i) the normal and tangential velocity 
boundary conditions are accurately satisfied, (ii) The zero divergence con- 
dition is accurately satisfied. On the other hand, in the original projection 
method, the step advancing the flow velocity forward in time enforces the 
velocity boundary conditions, but cannot guarantee that incompressibility is 
maintained. In the other step, the pressure is used to recover incompressibil- 
ity — by projecting the flow velocity field onto the space of divergence free 
fields. Unfortunately, while removing the divergence, this second step does 
not necessarily preserve the correct velocity boundary values il2|. 

Finally the modified equation ( 130|) implemented here differs from the pro- 
jection's method velocity forward step primarily by the boundary conditions 
imposed]^ Specifically, the projection method imposes Dirichlet boundary 
conditions for each component of the velocity field when propagating the 
flow velocity. Dirichlet boundary conditions have the obvious advantage of 



®In the projection method the contribution from the pressure is ignored in the momen- 
tum equation step, but we will not focus on this here — the prior two paragraphs deal 
with the consequences of this difference. 
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decoupling the velocity components. Numerically this means that an implicit 
treatment of the stiff viscosity operator dt — fJ^A is straightforward. By con- 
trast, on a regular grid, the boundary conditions in f pOj) couple the boundary 
velocities. Although this coupling does not pose a serious difficulty for ex- 
plicit schemes, it adds a (potentially serious) difficulty to the implementation 
of any scheme that treats the viscosity operator in the equations implicitly. 

6. Implementation 

The objective of this section is to study the convergence and accuracy of 
the scheme proposed in § |5l In particular: to verify the theoretical prediction 
that the scheme is second order in space, all the way up to the boundary. To 
this end, here we present the results from numerical computations done with 
two examples of flows driven by an external forcing f on a domain: (i) Flow 
on a square domain in § 16. ![ and (ii) Flow on an irregular domain in § 16.21 
In these examples the external forcing is selected so that an exact (analytic) 
solution to the equations can be produced In each case we compare 

the numerically integrated fields, u and p, with the known exact fields, and 
the errors are computed in an appropriate norm (see below). Note that, as 
pointed out at the beginning of § |5l here we only solve the linear Navier- 
Stokes equations. At some level, the nonlinear terms can be thought of as a 
non-zero divergence external forcing function. Consequently, to ensure a fair 
test, we choose f so that it is not divergence free. 

In what follows the numerical errors are measured using the discrete 
grid norm defined by 

\\g\\c^ = max{\gij\}. (44) 

Here the maximum is over all the indexes i and j — corresponding to the 
appropriate grid coordinates (xj, yj) in the computational domain — for 
the field g in the staggered grid. Namely: nodes for the pressure p, mid- 
points of the horizontal edges for u, and mid-points of the vertical edges 
for V. Note that neither ghost pressure points, nor boundary velocities, are 
included within the index set in f l44|) — we consider these as auxiliary nu- 
merical variables, introduced for the purpose of implementing the boundary 



^That is, proceed backwards: first write an incompressible flow, and then compute the 
forcing, and boundary conditions, it corresponds to. 



28 



conditions, but not part of the actual solution. Finally, for the exact (ana- 
lytic) solution, grid values are obtained by evaluation at the points (xj, yj). 
For example, to compute the error in the pressure, first define the discrete 
field Cij = Pij — p{xi, Uj) — where p is the numerical pressure field and p is 
the exact continuous field, and then compute the norm above for e. 

6.1. Flow on a square domain 

Here we present the results of solving the linear Navier-Stokes equations 
on the unit square < x, ?/ < 1, with no-slip and no flux boundary conditions 
(i.e.: u = t> = on dVL), and viscosity /i = 1. We selected the forcing 
function f = Uf + Vj5 — /iAuto yield the following (incompressible) velocity 
and pressure flelds: 



Using as initial conditions those provided by (l45l - H6|) at t = 0, we solved the 
equations for various grid sizes Ax, with time step At = 0.2 (Ax)^ — for 
stability: see § 15. equation (!32|) . Then we compared the numerical results 
with the exact flelds in f l45fH7|) . For instance, flgure [3] shows the error flelds 
for the velocity and pressure at time t = Air, for a (fairly coarse) 40 x 40 
grid. Note that the error is fairly uniform in size across the whole domain. 
Unlike the common situation with projection methods [l2[, there are neither 
numerical boundary layers, nor numerical corner layers. 

Since the algorithm is second order consistent in space, and flrst order in 
time, with At ~ (Ax)^, we expect the errors to scale with (Ax)^. This is 
precisely what we observe, with the errors measured in the L°° norm — see 
equation (jH]). We omit the convergence plots for the current case, and show 
them for the case presented in § 16.21 only. 

This example is not particularly taxing in terms of the algorithm imple- 
mentation, since the boundaries are parallel to the grid lines. Instead, the 
aim here is to illustrate the role of the feedback parameter A. For this pur- 
pose we present here the results of two tests. In the flrst test, we evolved the 
numerical solutions (with a flxed A x) for different values of A. To be precise, 
we used a 40 x 40 grid0 Figure m illustrates the results of this flrst test, with 



'With viscosity n — I and a a time step At = 0.2 (A x)'^. 



m(x, y, t) 
f (x, y, t) 
p(x, y, t) 



71 cos(t) sin(2 7r?/) sin^(7rx) 
TT cos(t) sin(2 7rx) sin^(7r?/) 
cos(t) cos(7rx) sm{ny). 



(45) 
(46) 
(47) 
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Figure 3: Error fields for tlie numerical solution corresponding to the exact 
formulas in f l45] - H7j) . The plots are at time t = Air, for a 40 x 40 grid, with 
At = 0.2 (A x)^. The horizontal velocity u (top) and the pressure p (bottom) 
error fields are shown. The error is uniform in size across the domain. There 
are neither numerical boundary layers, nor numerical corner layers. 

plots of the error in the horizontal velocity u, as a function of time, for 
various values of A. The errors in the pressure p, and in the vertical velocity 
V, exhibit the same qualitative behavior. In the absence of feedback — i.e. 
A = 0, the errors grow steadily in time, as can be expected from equation 
f ITS]) when influenced by the presence of numerical noise. Even though this 
effect does not constitute an instability (in the sense of exponential growth), 
there appears to be no bound to the growth. Thus, after a sufficiently long 
time, the errors can become substantial. The figure shows also that moderate 
values of A — i.e. A ~ 10 or bigger, are enough to control the errors. 

Physically, the errors that occur when A is too small correspond to fiuid 
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Figure 4: Time evolution of the L°° error in u, for the flow in a square 
domain, on a 40 x 40 mesh, with different values of the parameter A. For 
A = the error grows in time. Values larger than A ~ 10 control the error. 



leakage through the domain walls. However, these errors cause no net mass 
loss or gain, since V ■ u = applies — actually, V • u = O ((Ax)^), as we 
show later: see figure [TH Any positive flow across some part of the boundary 
must be compensated by a negative flow elsewhere. 

In the second test, we introduce artificial random errors to the normal ve- 
locity along the boundary. The purpose of this second test is to study, under 
a controlled situation, the error sensitivity of the normal velocity bound- 
ary condition implementation — the purpose for which the parameter A was 
introduced in § HI Specifically, the main concern are the errors that are in- 
troduced by the implementation of the Neumann boundary condition for the 
pressure on non-conforming boundaries (as happens for the example in § l6.2p . 
For simple geometries, such as a square, a finite differences discretization of 
the equation can be easily designed so that the discrete analog of the solv- 
ability condition applies. For curved geometries, however, this is generally 
not possible — therefore, a least squares solution of the equation is required, 
which introduces errors throughout the pressure field. Further, even though 
these errors are small (of the same order as the approximation scheme used) , 
they are hard to quantify in detail — e.g. write a leading order approxima- 
tion in terms of derivatives of the exact solution. This motivates this test. 



31 



which shows that the normal boundary condition calculation is robust. 



0.06 




t 



Figure 5: Time evolution of the error in u, for the flow in a square domain, 
on a 40 X 40 mesh, with a random forcing on the boundary. Without feedback 
(A = 0, solid line) the error grows steadily. With feedback (A = 100, dashed 
line) the error growth saturates at a small value. 

To perform the test, we introduce independent random errors at each 
point on the domain boundary, at the start of each Euler time step. To be 
precise, the boundary condition for equation ( 1311) is taken as 

n ■ Vp" = n ■ (f" + Au" + A u") + r, (48) 

where, at each boundary point r is sampled randomly from the interval [0, 1]. 
This represents an 0(1) random perturbation to the normal component of the 
desired velocity, g = 0, at the boundary]^ Hence, it is a very demanding test 
of how insensitive the boundary condition implementation is to errors. On 
the other hand, because of the highly "oscillatory" nature of the perturbation, 
the effect it has is to merely produce a thin boundary layer on the pressure 
— with the perturbations exponentially damped away from the boundary. 
Hence the perturbation only affects the solution at the boundary, and a large 
enough value of A can keep the solution under control — as we show next. 

Figure |5] shows the error in the horizontal velocity m as a function of time, 
both for the feedback controlled solution with A = 100, and the undamped 



Recall that in this section we take n — I, and neglect the nonlinear terms. 
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solution with A = 0. The errors saturate in the feedback solution, but they 
grow steadily for the undamped solution. Without damping, the drift errors 
contribute sizable effects to the solution after a short time period. 

Finally, figure [6] illustrates the flow leakage produced by the random er- 
rors. This figure shows a cross section of the horizontal velocity u after 1000 
time steps, for two values of the control parameter: A = and A = 100. The 
first value shows a significant flow through the boundary, while the second 
does not. This plot also illustrates the point made earlier: there is no signif- 
icant mass loss (gain) even when A = 0, with positive flows compensated by 
negative flows elsewhere. 
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Figure 6: Velocity cross-section u{x, y = 0.4872), at time t = 0.1315 (1000 
time steps), for the flow in a square domain, on a 40 x 40 mesh, with a 
random forcing on the boundary. The dashed line (A = 100) reproduces the 
correct field with zero flux at the boundaries. The solid line (A = 0) shows a 
non-zero flow through the boundary. 



6.2. Flow on an Irregular Domain 

In this section we implement the numerical scheme described in § [5l for 
an example with an externally forced velocity field on an irregular shaped 
domain. An exact solution in the domain, needed to compute the errors, is 
constructed in a similar fashion to that in § 16.11 However, instead of the unit 
square, the selected domain Q is the 2x2 square < x, y < 2, with the 
disk of radius 1/4 centered at (3/4, 1) removed. In addition, we use periodic 
boundary conditions (with period 2) in the y direction, and impose a nonzero 
flux (u = g 7^ 0) on the rest of the boundary. 
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To construct an incompressible velocity field in the domain f2, we use 
a stream function ip, with u = {u, v) = {ipy, —ipx)- Then we prescribe 
a pressure, and choose the forcing function by f = uj + Vp — /xAu, with 
viscosity = 1. The function g, used to prescribe the boundary condition 
for the velocity, is selected to match the values given by {ipy, —ipx)- 

In this example we set the feedback parameter A to A = 100, considerably 
larger than the minimum needed to get good behavior in the calculations in 
§ 16.11 This is to be expected from the results of the last test in § 16.11 — 
see figures [5] and [61 Because of the curved inner boundary in the current 
example, errors in the implementation of the Neumann boundary condition 
for the pressure occur, which trigger a drift in the normal boundary condition 
for the velocity — unless a large enough A is used. Thus, for example, values 
where A = 0(10) did not allow the feedback to quickly track fluctuations in 
the normal velocity component of the boundary conditions. 



Figure 7: Contour lines for the stream function ip, defined in equation (|50|) . 
on the irregular domain of the second example. The domain is a 2 x 2 square, 
with a hole of radius 1/4 centered at (x, y) = (3/4, 1). The stream function 
is periodic, of period 2, in the y-direction. 

The stream function ip = ipi^^ Z/? t) (see figure [7]) is defined as follows: 
First, introduce the function i/jq = ipoix, y), defined on R^, by 




2 



X 




(49) 
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where r = ^/{x — 3/4) + {y — 1)^. Then construct a 2-periodic function in y 
by adding shifted copies ipo, as follows 



oo 



ip{x, y, t) = cos 



(50) 



Finally, the pressure is given by the formula 



p{x, y, t) = sin(t) cos(7rx) sin(7ry) 



(51) 



Figure [7] shows the contour lines for on the domain Q. Note that, although 
both ipo and ^ ^/jq vanish at the edge of the circle centered at (3/4, 1) — the 
inner boundary of Q, the shifted copies of i/jq do not. Therefore, the resulting 
flow velocity [ipy, —ipx) has some small components, both in the normal and 
tangent directions to the circle. 



To illustrate the second order convergence of the scheme, below we present 
the results of evolving the velocity and pressure flelds (for the problem de- 
scribed above), for varying grid sizes Ax, up to the flxed time t = 0.0657, 
using A = 100 and At = 0.2 (A x)^. Note that the selected flnal time ap- 
pears small, but it is large enough to require a number of time steps in the 
range O(IO^) to 0(10"^) — for the grid sizes we used. This is large enough 
to provide a reliable measure of the order of the scheme. 

Figure M shows the convergence of velocity and pressure, while flgure [9] 
shows the convergence of the partial derivatives of the horizontal velocity 
u. The flgures indicate second order L°° convergence for the velocity u, the 
pressure p, and even the velocity gradient {u^, Uy) — see § 16.31 



We illustrate the (typical) error behavior — both in time and in space, by 
showing the results of a calculation done at a flxed resolution. Speciflcally, 
we take an 80 x 80 grid, with At = 0.2(Ax)^ and A = 100, and solve the 
forced system of equations in this "flow on an irregular domain" example — 
recall that yU = 1. 

Figure [10] shows the time evolution of the L°° (spatial) errors in the 
pressure and in the horizontal velocity. It should be clear that, while the 



Second order convergence. 



Typical error behavior. 
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Figure 8: (Flow on an irregular domain example). Convergence plot for the 
pressure (circles: o) and the horizontal velocity (squares: □). The errors 
(in the L°° norm) are computed at the fixed time t = 0.0657, for different 
grid resolutions, with At = 0.2 (Ax)^. The slope of the solid straight line 
corresponds to the second order scaling error oc (Ax)^. 

errors oscillate (over almost one decade in amplitude), they do not exhibit 
any measurable growth with time. 

Figure [11] shows the horizontal velocity field at time t = 4.25 vr, while 
figure [12] shows the associated error field. Notice that, while the error is 
largest near the internal boundary — as expected from the difficulties that a 
curved, non-conforming, boundary causes — it is fairly well behaved, without 
abrupt transitions on the scale of the grid size. This explains why the error 
in the gradient of the velocity is also second order accurate — see figure [9] 
and remark [9] a feature that should be particularly useful for the calculation 
of fiuid-solid forces (stresses) along domain boundaries. 

Similarly, figure [12] shows plots of the pressure and of the associated 
pressure error field, again for the time t = 4.25 tt. Just as in the case of the 
velocity, the errors are dominated by what happens at the internal circular 
boundary. The error near the internal wall is a bit more jagged than the 
error for the velocity field. We did not check the order of convergence for the 
gradient of the pressure since one does not need the gradient of the pressure 
to compute forces. 

Finally, figure [H] shows a plot of the error in the numerical divergence of 
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Ax 

Figure 9: (Flow on an irregular domain example). Convergence plots for the 
partial derivatives of the horizontal velocity: Uy (circles: o) and (squares: 
□). The errors (in the L°° norm) are computed at the fixed time t = 0.0657, 
for different grid resolutions, with At = 0.2 (Aa;)^. The slope of the solid 
straight line corresponds to the second order scaling error oc {Axf. The 
fact that the errors for the gradient of the velocity are also second order is 
important — see § 16. 3[ 



10-^ 




^° 2 4 6 8 10 12 14 

t 

Figure 10: (Flow on an irregular domain example, on an 80 x 80 grid). 
Evolution in time of the L°° norms of the (spatial) errors for the pressure 
(dashed curve) and horizontal velocity (solid curve). 
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Figure 11: (Flow on an irregular domain example, on an 80 x 80 grid). 
Numerical horizontal velocity field u at t = 4.25 vr. 



X 10"" 




Figure 12: (Flow on an irregular domain example, on an 80 x 80 grid). Error 
field for the horizontal velocity u at t = 4.25 vr. 

the flow field, also for time t = 4.25 tt. The errors in the divergence are also 
second order. However: notice that they are much smaller than the errors 
in the flow field or in the pressure. This is indicative of the rather strong 
enforcement of incompressibility that equations Q) and (fTO!) guarantee. 
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Figure 13: (Flow on an irregular domain example, on an 80 x 80 grid). Left: 
numerical pressure field p at t = 4.25 vr. Right: associated error field. 




Figure 14: (Flow on an irregular domain example, on an 80 x 80 grid). 
Numerical error in the divergence of the velocity at t = 4.25 tt. The O(10~^) 
amplitude is the result of the 0{Ax^) order of the scheme. 

6. 3. Convergence of the derivatives 

A key reason to worry about uniform convergence, up to the boundary, of 
the numerical solution is that this is a necessary condition for the accurate 
modeling of solid-fluid interactions. The evaluation of fluid-solid stresses 
requires both the pressure, as well as the derivatives of the velocities, to be 
accurate at the boundaries. Hence, in this paper we investigated the behavior 
of the errors for not just the pressure p and the velocity fields (u, f ), but for 
the gradient of the flow velocity as well. 
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An important point to notice is that figure M shows that 

The derivatives of the velocities exhibit errors that 
appear to be second order as well — in the norm. 

At first sight this may seem surprising: since the velocities are second order 
accurate, one would expect their derivatives to be first order only. However, 
this is a worse case scenario, based on the assumption of random errors. 
On the other hand, the errors for finite differences approximations (on reg- 
ular grids) are typically not random: If the solutions are smooth enough, 
the errors can be expanded in powers of Ax, with coefficients that involve 
derivatives of the solution. Hence, in this case, taking low order derivative^ 
does not degrade the order of convergence. This, of course, is an important 
advantage of finite differences over other approaches (e.g. projection) which 
do not have this property. 

The astute reader may have noticed that the argument in the prior para- 
graph avoids the issue of boundary condition implementation, which can ruin 
the smooth error expansions available for finite differences. This requires 
some extra considerations: 

^ On conforming boundaries, it is usually possible to approximate the 
boundary conditions in such a way that smooth expansions (in powers 
of A x) remain available for the truncation errors. 

On non-conforming boundaries, on the other hand, it is very easy to 
ruin the smooth error expansions. The local stencils used to approx- 
imate the boundary conditions along the boundary will, typically, ex- 
perience abrupt changes in response to the placement of the boundary 
relative to the regular grid. 

This is where the global approach that we used in § 5.2 to implement the 
boundary conditions comes to the rescue: by linking each local stencil to its 
neighbors, the non-smoothness caused by abrupt stencil changes is smeared. 
Thus a better behaved error is obtained, which then explains why the result 
in flS2]) occurs. Notice that, as of now, we neither know if the smearing 
process is enough to make the errors at leading order — which is what 



^^At high order, round off errors dominate. There is an irreducible 0{e {Ax) contri- 
bution from them, where e is the round off error, and q is the order of the derivative. 
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is needed to make the errors in the velocity gradient second order, nor if 
the errors in the velocity gradient are actually second order. However, the 
numerical evidence seems to point in this direction. 

A final point is that, the "simplest" way to get around the issue in item [6b 
above, is to implement the boundary conditions at a higher order than re- 
quired. This has the disadvantage that it can lead to a messier than needed 
algorithm, but (generally) it should not significantly increase the computa- 
tional cost — since it only involves the boundary points of the grid. 

7. Conclusions 

Through the introduction of the pressure Poisson equation with consistent 
boundary conditions, we give an equivalent formulation of the incompress- 
ible Navier-Stokes equations. In this formulation the momentum equation 
takes the form of a vector heat equation with unconventional boundary con- 
ditions, while the pressure Poisson equation can be used to explicitly describe 
the pressure as a function of velocity at any fixed time. It follows that the re- 
formulated system of equations is ideal for using efficient numerical marching 
methods, where there are no particular theoretical limitations to the order 
of accuracy or the method of implementation. 

In addition, we devise and implement a second order discretization (uni- 
form up to the boundary) of the equations on irregular domains. We address 
the issue of numerical stability for the normal boundary velocity by adding 
an appropriate feedback term to the equations. For the scheme, we use finite 
differences on a regular, staggered grid, so that the domain boundary can be 
immersed within an x M mesh. We discretize to second order spatial ac- 
curacy, and verify the order of accuracy in L°°, both for the velocity and the 
pressure on an irregular domain. Hence, the solution^ converge 0{{AxY) 
uniformly (all the way up to the boundary) as the grid spacing Ax — )■ 0. 
Although, we formulate and implement the scheme in two dimensions, the 
algorithm has a natural extension to three dimensions. 

There are several issues and extensions that we hope to address in future 
work. First, for irregular domains on a regular grid, the proposed momen- 
tum equation boundary conditions implicitly couple all the components of 
the velocity field. Therefore, the resulting vector heat operator, dt — fJ,A, 



■^Note that we have explored only the case where the solution is at least C*. 
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cannot be solved "component by component". This makes the implementa- 
tion of schemes that do a naive implicit treatment of the stiff viscosity term 
computationally expensive, is there a better way? Second, since we handle 
the pressure boundary conditions implicitly, our discretized Poisson equa- 
tion does not have a symmetric matrix — even though the deviations from 
symmetry are "small" (i.e.: just the rows involving the boundary terms). In 
practice, the Poisson inversion dominates the computational expense. Hence 
a symmetric formulation would be very desirable, both for numerical sta- 
bility reasons, as well as to be able to exploit fast Poisson matrix solvers. 
Third, the formulation in this paper is for fixed domains, but extensions to 
deformable and or moving geometries seem possible. Fourth, can the ideas in 
this paper be extended to flows with variable densities (e.g.: stratified flows)? 
Lastly, there is the issue of obtaining implementation of the boundary con- 
ditions that yield smooth errors. As explained in § 16. 3[ this guarantees that 
the derivatives of the flow velocity have the same order of accuracy as the 
velocity itself, allowing an accurate calculation of the fluid forces on objects 
immersed in the flow (as well as the domain boundaries). 
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Appendix A. Further Modification for Solvability 

In this appendix we address the question posed and motivated in re- 
mark [71 Namely: Can the equations in l[2U^21\) be modified in such a way 
that they make sense even for initial conditions that are not incompress- 
ible? In fact, in such a way that if a solution starts with V ■ u 7^ and 
n ■ (u — g) 7^ 0, then (as t — )■ 00) V ■ u — )■ and n ■ (u — g) — )■ — so that 
the solution converges towards a solution of the Navier-Stokes equation. 
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As pointed out in remark [71 the problem with (12UH2T]) is that the solv- 
abihty condition for the Poisson equation (121 p is not satisfied when V-u 7^ 0. 
Hence the equations become ill-posed, as they have no solution. The obvious 
answer to this dilemma is to interpret the solution to fl2Tl) in an appropriate 
least squares sense, which is equivalent to modifying the non- homogeneous 
terms in the equation by projecting them onto the space of right hand sides 
for which the Poisson equation has a solution. Symbolically, write fl^Tj) in 
the form 

Ap = g for X G I 
n ■ Vp = h for x G d^l, J 

where g and h are defined in (12T|) . Then modify the equation to 

= 9p fo^^ X ^ ^) 
n ■ Vp = hp for x G d^l, 

where {gp, hp) = P (^f, h) for some projection operator P such that 

[ gpdV= [ hp dA. (A.3) 
Jn Jan 

The question, however, is: which projection? 

The discretization failures in solvability arising during the course of a nu- 
merical calculation are small, since it should be {ge, he) — {g, h) = O ((Ax)'') 

— where q is the order of the method, and {ge, he) is the exact right hand 
side. Thus one can argue that, as long as {gp, hp) — {g, h) = 0((Ax)''), the 
resulting numerical solution will be accurate to within the appropriate order 

— see remarks [7] and [TOl In this section, however, the aim is to consider 
situations were there is no small parameter (i.e. A x) guaranteeing that solv- 
ability is "almost" satisfied. In particular, we want to consider situations 
where V-u ^ 0, and |V ■ u| ^ 1 does not apply — leading to errors in 
solvability which are not small. It follows that here we must be careful with 
the choice of the projection. 

Obviously, a very desirable property of the selected projection is that it 
should preserve the validity of equations fP- nU]) — so that the time evolution 
drives V ■ u to zero. Hence: it must be that gp = g, with only h affected 
by the projection. Further: since the solvability condition involves h only 
via its mean value over dQ, the simplest projection that works is one that 
appropriately adjusts the mean of h, and nothing else. Thus we propose to 
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modify the equations in (12D] - [^T]) as follows: leave f l20p as is, as well as the 
imposed boundary condition constraint dl]), but replace ( 1?T]) by 

Ap = -V- ((u- V)u) + V-f for xG Q, 

n ■ Vp = n ■ (f - gt + A u - (u ■ V) u) } (A.4) 

+ A n ■ (u - g) - C for xedn, 



where 



^ / n ■ Au + Au) dA, (A.5) 
^ Jan 



and S = Jq^ d A is the surface area of the boundary. In terms of flA.lHA.2|) 
this corresponds to the projection 

gp = g and hp = h — C, (A. 6) 

where 

C = ^ (^J^^hdA- JjdVy (A.7) 

This is clearly a projection, since C = for {gp, hp), so that = P. Further, 
since the solvability condition for ( \A.l\i is precisely C = — see equation 



(lA.Sp . the solvability condition for ( \AA\i is satisfied — even if V • u 7^ 0, 



though (of course) C = when V ■ u = 0. 

Finally, we remark (again) that the projection in (1A.6P is not unique. In 
particular, numerical implementations of the Poisson equation with Neumann 
conditions often use least squares projections, which alter both the boundary 
condition h and the source term g. This makes sense if the solvability errors 
are small. However, in general it seems desirable to not alter the source term, 
and keep flU HTU]) valid. This still does not make (IA.6P unique, but it makes 
it the simplest projection. Others would also alter (in some appropriate 
eigenfunction representation) the zero mean components of h. This seems 
not only un-necessarily complicated, but it may also affect the validity of 
the result stated below equation (lA.Sp . defeating the whole purpose of the 
reformulation in this appendix. 

The system of equations in (1201) and (lA.4p makes sense for arbitrary 
flows u, which are neither restricted by the incompressibility condition (f) = 
V ■ u = in f2, nor the normal velocity boundary condition £^ = n ■ (u — g) = 
0. Furthermore: this system, at least in bounded domains Q, includes the 
Navier-Stokes equations as a global attractor for the smooth solutions. This 
is easy to see as follows: 
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First, because of equations (P-[TUI). decays exponentially, at a rate con- 
trolled by the smallest eigenvalue of L = —A in Q, with Dirichlet boundary 
conditions. In particular: 




an 



n-(;uAu + Au) dA 



Jn 



[ (/iA0 + A0) dV 



vanishes exponentially. 

Second, it is easy to see that, for the system in (120|) and (lA.4p . equa- 
tion (I22!) is modified to 



where £ = n • (u — g). Hence S also vanishes exponentially. 

Of course, if = and £ = initially, then they remain so for all times, 
and the evolution provided by (|2UI) and (lA.4p is, exactly, the Navier-Stokes 
evolution. 

In conclusion, the formulation in this section is not only an interesting 
theoretical fact. It also provides a robust framework within which numeri- 
cal solvers for the incompressible Navier Stokes equations can be developed, 
without having to worry about the (potentially deleterious) effects that dis- 
cretization (or initial condition) errors, can cause when they violate mass 
conservation — because either V ■ u = in fi, or n ■ (u — g) = in fail. 
In addition, the formulation eliminates the necessity of having to enforce the 
condition V ■ u = directly, which is a core difficulty for the solution of the 
Navier-Stokes equations. 

Appendix B. The V • u = Boundary Condition 

In the main body of the paper, the V ■ u = boundary condition is 
implemented using a regular Cartesian grid. In Cartesian coordinates, this 
boundary condition relates the horizontal and vertical velocities (for two 
dimensions) via the standard formulas 



= -X£ + C for xG 



(A.8) 



u 



u i + V j, 
du dv 
dx dy' 




Vu 



(B.2) 



where i and j are the coordinate unit vectors. 
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For some applications, it may be more convenient to implement the nu- 
merical computation on a coordinate system where the boundaries are con- 
forming. In general curvilinear coordinate systems, the divergence acquires a 
more complicated form than the one above. Hence, it is our purpose here to 
display the form that the V ■ u = boundary condition takes for a conform- 
ing boundary in a general orthogonal curvilinear coordinate system. We do 
the 2-D case only — the 3-D case is quite similar. 

Appendix B.l. Curvilinear coordinates and a conforming boundary 

In some region near the boundary dVt oi the domain of integration, as- 
sume that an orthogonal, curvilinear, set of coordinates (77, ^) has been se- 
lected — such that the boundary is given by 77 = c, for some constant c. In 
terms of the coordinates (77, ^), the vector field u can then be written in the 
form 

u = t) + I, (B.3) 

where and are the velocity components in the coordinate directions 
given by the unit vectors i) and ^, respectively. Let now the functions 
X = x{ri, C,) and y = y{ri, ^) describe the relationship of the curvilinear 



coordinates with a Cartesian system. Then (e.g. see 151 ]) 



1 f dx ^ dy \ . - I f dx ^ dy \ 
where and are local normalization factors 

Further, re-interpreting and as functions of (r/, ^), we have 

V-u = -i- + (B.6) 

The V ■ u = boundary condition is accompanied by the Dirichlet boundary 
condition n x u = n x g for the vector momentum equation. For a conforming 
boundary, this is just = along the boundary curve where rj = c. Thus 
the V ■ u = boundary condition reduces to a Robin boundary condition on 
the normal velocity (provided that the boundary is piecewise smooth) 

d d 

j^{s^u,,) + -^{s^g^) = (B.7) 

along rj = c. 
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Appendix B.2. Example: polar coordinates 

Consider the special case of a polar coordinate system x = rcos6 and 
y = rsin^^, with the boundary at r = 1. Let the angular velocity at the 
boundary be given by Ug{l, 9) = ge{d). Then since (in this case) Sg = r and 
Sr = 1, the zero divergence boundary condition becomes 

-^{rur) + ^9e^0 for r = 1. (B.8) 
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